Excitonic Instabilities and Insulating States in Bilayer Graphene 
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The competing ground states of bilayer graphene are studied by applying renormalization group 
techniques to a bilayer honeycomb lattice with nearest neighbor hopping. In the absence of interac- 
tions, the Fermi surface of this model at half-filling consists of two nodal points with momenta K, 
K', where the conduction band and valence band touch each other, yielding a semi-metal. Since near 
these two points the energy dispersion is quadratic with perfect particle-hole symmetry, excitonic in- 
stabilities are inevitable if inter-band interactions are present. Using a perturbative renormalization 
group analysis up to the one-loop level, we find different competing ordered ground states, includ- 
ing ferromagnetism, superconductivity, spin and charge density wave states with ordering vector 
Q = K — K', and excitonic insulator states. In addition, two states with valley symmetry breaking 
are found in the excitonic insulating and ferromagnetic phases. This analysis strongly suggests that 
the ground state of bilayer graphene should be gapped, and with the exception of superconductivity, 
all other possible ground states are insulating. 



I. INTRODUTION 

Graphene is a quasi-2D carbon material with 
a honeycomb lattice structure. Its band struc- 
ture is captured by a tight binding model, as 
illustrated in Fig. [TJ with two interpenetrating 
triangular sublattices a and b 
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where denotes a sum over all nearest neigh- 
bor pairs. At the charge neutrality point, this 
model yields a semi-metal for which the Fermi 
surface (FS) contains only two nodal points. 
Since the energy dispersion is linear in the vicin- 
ity of these Dirac points, the corresponding low- 
energy effective Hamiltonian is given by a 2D 
Dirac model. This unique electronic structure 
leads to many interesting phenomenaji 

Although interactions between electrons are 
present in graphene, the one-particle picture 
works surprisingly well. In contrast to ordi- 
nary metals, the ground state of the electrons 
in graphene does not behave like a Landau 
Fermi liquid, but rather belongs to the univer- 
sality class of Dirac liquids.— One of the differ- 
ences between these ground states is that short- 
range interactions between electrons are irrele- 
vant in Dirac liquids £ This may explain why 
the one-particle picture is applicable, regard- 




FIG. 1. Bilayer graphene with AB-stacking: oi, i>i 
are the two sublattice sites in the upper layer, a-2, 
62 are the two sublattice sites in the lower layer. 
70 is the tight-binding hopping constant between 
01; 61, 71 is the hopping between a\ and a% 73 is 
the hopping between b\ and 62- ai = |(3, \/3) and 
&2 = f (3, — \/3) are the primitive lattice vectors. 



less of the perfect particle-hole nesting prop- 
erties of the lattice. However, recent experi- 
ments have shown evidence that the Dirac cone 
is renormalizcd,— suggesting that electron inter- 
actions are important on some level. Recently, 
the interactions between electrons in graphene 
have been modeled by a long-range Coulomb in- 
teraction or by using an effective (2+l)D QED 
modelj£££ 

For bilayer graphene (BLG), tight-binding 
calculations also show that the non-interacting 
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ground state is a semi-metal. But in this 
case, the dispersion near the FS points is 
quadratic rather than linear^ Because of this, 
all short-range interactions now become rele- 
vant perturbations, and recent theories have 
predicted various possible spontaneous symme- 
try breaking ground states^ - — Furthermore, re- 
cent experiments^— have shown some evi- 
dence for FS reconstruction in BLG. These find- 
ings contradict the simple one-particle picture 
for BLG, based on a tight-binding model, and 
rather suggest that interactions between elec- 
trons play an important role in breaking down 
the FS. 

In this paper, the instabilities in BLG will 
be addressed by using a perturbative renormal- 
ization group approach. We consider the bi- 
layer honeycomb structure with nearest neigh- 
bor hopping as the low-energy effective model 
for BLG. Particle-hole symmetry is assumed, 
and RG arguments are used to identify the 
dominant channels and eliminate the irrelevant 
channels due to the interactions in the model. 
Using this setup, an array of possible ordered 
phases is found, which are competing with each 
other. In the following sections, the details of 
the model and the results and implications of 
our calculations will be discussed. 



II. BILAYER GRAPHENE AND THE 
MODEL HAMILTONIAN 



The crystal structure of BLG is given by a 
Bernal AB stacking of two sheets of graphene, 
shown in Fig. [TJ. In the absence of interactions, 
its band structure is effectively described by a 
tight-binding modeli In momentum space, the 
one-particle Hamiltonian with 74 ~ is given 
by 




(a) First Brillouin Zone (b) Energy dispersion 
around the Fermi 
points 



FIG. 2. (a) K = §1(1,-1=) andK' = 
are the two points, constituting the Fermi surface 
of the non-interacting system. A < 71 is the energy 
cutoff of the theory, dA is a thin shell contain high 
energy modes to be integrated out. (b) e c (K) and 
£f(K) are the dispersion energy of the conduction 
and the valence band respectively. The other two 
bands are gapped by 71 . 
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= (^^U^L^L). f( K ) = 

e %K-Si are ^ e orbital field operators, and 

f(l,V3), = f(l, -V3), S 3 = o(-l,0) 
are nearest-neighbor in-plane displacement vec- 



tors (a is the lattice constant). Fig. |2(a) shows 



the 1st Brillouin zone in momentum space with 
reciprocal vectors 61 = |^(l,\/3) and 62 = 

£(1 -V§). 

Since only low energy excitations are of inter- 
est here, we expand f(K) near K and K' (up 
to a phase factor e OT / 6 ), 



f(K) 



3 a 



k at K, f{K] 



3 a 



k* at K', 



where k = k x + iky, k = (k x ,k y ) is a small 
momentum deviation from K, K', and |fc| < 
A<|K|,|K'|._ 

In the following discussion, the trigonal warp- 
ing term 73 will be neglected. (The justification 
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for this will be discussed in the Sec. IVI[) . As ten as 



shown in Fig. 2(b) the resulting tight-binding 
band structure consists of 4 bands. Two of these 
bands arc gapped by 71 from the FS, whereas 
the other two bands touch each other at the K 
and K' Fermi points. This is similar to single- 
layer graphene, but for the bi-layer case the en- 
ergy dispersion is quadratic at the Fermi sur- 
face, 



e cJ (K)~±^-k 2 atK,K'. 



7i 



(2) 



In the following analysis of instabilities, the 
gapped bands will be ignored, because they are 
not important in the low energy limit. Be- 
fore writing down the model Hamiltonian, let 
us introduce the creation (annihilation) opera- 
tors for electrons in bands e c {K) and tf(K) to 

be c K a (c\ <a ) and f Kc (f Ko .) respectively. c Ka 
and ]k<j are linear combinations of the local or- 
bital field operators (bi Ka , a 1K a, a 2 Kcr 7 b 2K<7 ): 

CKa = C C K ■ $ Ka , f Ka = C f K ■ V Ka , (3) 



where C K 



(Ci K ,Ci K ,Ci K ,Cl K ). These co- 
efficients near the Fermi surface can be found in 
Ref. Note that if the model is written us- 
ing the local orbital basis, in momentum space 
these coefficient C^- account for the form fac- 
tors in the interaction terms of the Hamiltonian. 
Because of the small fc dependence in the form 
factors, this can complicate the RG analysis. In 
order to avoid this problem, it is natural using 
the Bloch wave basis to build an effective Hamil- 
tonian of BLG. Then, the non- interacting part 
of the model Hamiltonian can be represented as 



Ho 



K 



e c (K)cl K c aK + e f (K)fl K U K , (4) 



where summing over all a is implicitly assumed. 

Turning to the interaction part of the Hamil- 
tonian, we follow the approach outlined in 
Refl23l (for more details see Appendix |A"| . and 
require particle-hole symmetry of exchanging 
the valence and conduction bands. Then the 
electron-electron interaction term can be writ- 



flint - o E 

K 3 ,Ki 

{ U (K 3 K A K 2 Ki ) C^ 3(r a ,CK 2 a'C Kl a 

+ U 2 (KaKiKzK^c^J^, f K2 <j>c Kl * 
+ U 3 (K 3 K 4 K 2 K 1 )fl^ Kia ,f K , a ,c Kia } 
+ {exchange (c <-> /)}, 

where momentum conservation is implicitly 
contained in U (see Appendix |A|) . Here, the cou- 
pling constant Uq denotes the intra-band inter- 
action, whereas U\, IA 2 and U3 are inter-band 
interactions. 

So far, no explicit advantages are obvious 
by using the Bloch wave basis. In addition, 
the momentum dependence in the coupling con- 
stants complicates the study too. However, this 
complication will be removed due to the trivial 
topology of the FS. (See below) 



III. RENORMALIZATION GROUP 
ANALYSIS OF THE BLG MODEL 

Here, we apply the pertubative renormaliza- 
tion group (RG) method to explore the low- 
energy physics of the BLG model in the pres- 
ence of interactions, following the standard pro- 
cedure outlined in Ref. 

From a tree level analysis (see Appendix |B|) , 
we find that only a finite set of coupling con- 
stants are marginal. In the low energy limit, 
only the interacting channels which depend on 
K,K' are not rcnormalized to zero. The corre- 
sponding bare coupling constants arc listed and 
classified into Table HI 

Here, the subscripts 0, 1, 2 of the coupling 
constants indicate the various scattering pro- 
cesses between valleys. The difference between 
processes with 0, 1 versus 2 is that after scat- 
tering processes with 0, 1 do not exchange val- 
ley indices between two particles, but processes 
with 2 do. Therefore, the scattering processes 



4 



Uo lh U 2 U z 

W(K,K,K,K) h g u v 
W(K,K',K',K) hi gi ui vi 
U(K',K,K',K)\h 2 gg m 2 «g 

TABLE I. Bare coupling parameters of the 
marginally relevant processes. 



with subscript 2 always involve large momen- 
tum transfers. 

Since these coupling constants are marginal, 
performing one-loop corrections to the RG flow 
equations is necessary. Since the interaction is 
quartic, i.e. involving only two-body scattering, 
there are only three distinct channels to transfer 
momentum. Following the terminology of Ref. 
[H, these processes are named ZS, ZS', and BCS. 




FIG. 3. Feynman Diagrams: 1, 2, 3, 4 represent the 
low energy modes with momentum if 1,2,3,4, band 
index I, valley index a, and spin a. The momentum 
inside the loop, K, must lie within the shell dA, and 
Q = K 3 - Ki, Q' = K A - Ki, P = K x + K 2 . Note 
that the interaction lines are suppressed. 



The corresponding Feynman diagrams are 
schematically shown in Fig. [3] All modes in 
the loop are high energy and need to be inte- 
grated out. After rescaling back to the original 
phase space volume, the coupling constants are 
modified, i.e. they are flowing in a 12 dimen- 
sional space of couplings. 

In order to have non-vanishing one-loop cor- 
rections, in the ZS and ZS' diagrams the two 
propagators in the loop must pair up with a dif- 
ferent band. For BCS, both propagators must 
pair up within the same band. Those graphs 
that do not satisfy the above criteria contain 
double poles in the frequency uj contour inte- 
gration. With this, many contributions of these 
diagram can be eliminated, thus greatly simpli- 
fying the calculation. 

In this work, we consider flow equations for 
the couplings up to the one-loop level. Cumu- 
lant expansion and Wick's contraction are used 
in the calculation. This method is convenient 
to keep track of the prefactor for each different 
diagram. 

The loop momentum integration (bubble di- 
agram) can be evaluated, 

r 27T r A dOkdk 1 dt 
Jo JA- d A (2tt) 2 2| C/ (J0I " W 

where v§ = and dt = [ s the RG run- 

ning parameter. Therefore, the RG flow rate 
equations under one-loop correction are given 
by, 
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If go = gi = g-2 = 0, these RG flow rate equation can be solved exactly, and decoupled 

into a simple result, 
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((u Q - 2v ) - (u 2 - 2vi)) 2 
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(8) 



Before finishing this section, we need to ad- 
dress the effects of quadratic perturbations. 
The two most relevant perturbations are the 
chemical potcnal and trigonal warping, i.e. the 
73 hopping term. These perturbations are in 
principle relevant under the tree level, i.e. scal- 
ing as s 2 and s respectively. The chemical po- 
tential determines the density of the system, 
and the trigonal warping splits the original two 
Fermi points into four. 



However, the divergences of the susceptibil- 
ities (see next section) emerge at some finite 
energy scale, and the RG flow must be stopped 
at this point. This energy scale determines the 
ordered state mean field transition temperature 
T c . If T c is far above the trigonal warping recon- 
struction energy, this quadratic perturbation is 
not significant. This introduces an infrared cut- 
off to the validity of the analysis J^—. In addi- 
tion, the divergences also imply that the original 
FS is unstable towards opening a gap. Although 
one should follow the procedure in Ref. to 
fine-tune the chemical potential to keep the sys- 
tem density fixed, not carrying out this proce- 
dure does not affect the results significantly. 



Because of this, we argue that trigonal warp- 
ing and the chemical potential do not play an es- 
sential role in the analysis at the one-loop level, 
as long as the energy scale of the instabilities is 
found to be far beyond the infrared limit. 



IV. SUSCEPTIBILITIES AND 
POSSIBLE GROUND STATES 



To gain further understanding into the 
physics of BLG, we introduce test vertices into 
the original Hamiltonian i 25 i 26 These test ver- 
tices correspond to the pairing susceptibilities, 

E Q fc c L,fc T L' ® < s 'Jw,fc, (9) 



A SC" Eat 



+f as ,kT* a , <E>a v ss ,f, 



a' s' . — k 
a's' ,—k 



-fas,k.T* a , ®<J V ss ,l 



a' s' 



(10) 

(ii) 



where, \i,v = 0, x, y, z, t° = a are 2x2 identity 
matrices, t x ' v,z and <j x ' v ' z are 2x2 Pauli matri- 
ces, r denotes the valley degree of freedom with 
basis (K,K'), and a denotes the spin degree of 
freedom, j indicates the different pairings listed 
in Table M 

Performing an RG analysis at the one-loop 
level with this additional new perturbed Hamil- 
tonian, the vertices (As) are renormalized, and 
the new renormalized vertices are of the form 



1 



Af" = A,(I + i --r,ln S ), 
where the T,- are listed in Table HI 
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j 


Ordered State 




r i 


FM 


ferromagnetism 


r°0a z 


+ go + U2 + g2 


FM' 


FM without valley symmetry 


T z ® CT Z 


MO + go - («2 +fl2) 


SDW 


spin density wave 


r x ®a z 


"I + gi 


EI 


excitonic insulator 


r°®a° 


(uq - 2v ) + (U2 - 2«i) + (go - (52 - 2gi)) 


EI' 


EI without valley symmetry 


t z ®a° 


(uo - 2« ) - (u2 - 2«i) + (30 + {gi - 2ffi)) 


CDW 


charge density wave 


t x ®a° 


mi + ffi — 2(u 2 + g%) 


SC 


superconductor 


t x <g> a v 


-{hi -h 2 + {gi - gi)) 


sc 


superconductor 


T x ® a v 


-{hi -h 2 - {gi - gi)) 



TABLE II. LPairing susceptibilities corresponding to the competing ground states in the presence of inter- 
actions. 



A. Case I: go — gi — g 2 = 

If go = g% = .92 = 0, Eq. ((SJ) decouples the 
susceptibilities, and one obtains 



rf=°(t) 



r g=0 



(0) 



1 - j^-rr °(o)t 



(13) 



Whether and where the susceptibilities 
{Tj In s, where t = In s) diverge is determined by 
the bare coupling constants. Each divergence 
in T 9 j~°{t) indicates that the system has a ten- 
dency toward the corresponding ordered state, 
labeled by 'j'. The first instability in a given 
channel represents the most dominant ordered 
state of the system at low energy. 

For the case go = g\ = 32 = 0, the situa- 
tion is relatively simple. If only repulsive in- 
teractions are considered, FM and SDW are 
the dominant instabilities. In order to produce 
instabilities in the other channels, fine-tuning 
of the bare parameters is needed. The 1^(0) 
must be positive such that other mean field so- 
lutions exist. More generally, since the parame- 
ter space spanned by the bare couplings is very 
large, constraining the search is desirable in or- 
der to make the exploration and analysis of the 
phase diagram meaningful. 

The relative strength between the bare cou- 
plings can be estimated. In general, the scat- 
tering processes within the same valley ho, go, 
uq, vq and hi, gi, ui, vi arc expected to be 
larger than h%, 52, U2, v%, because intra- valley 



scattering processes involve only small momen- 
tum transfer.— Applying these constraints, in 
Fig. |4] we show how these pairing susceptibili- 
ties compete with each other for a representa- 
tive choice of bare coupling parameters. In this 
example, the dominant low-energy divergence 
occurs in the FM channel, followed by FM' , 
SDW, CDW and SC at higher energy scales. 




FIG. 4. Flow of the susceptibilities: Here, we set 
ho = uo = vo, hi = in = vi = O.8/10, hi = 112 = 
v 2 = O.l/io, and go — gi — g 2 — {ho > 0) . The 
FM instabilities occupy a large region in parameter 
space, and fine-tuning is not necessary. 

The instabilities FM, FM', and SDW in- 
dicate broken spin symmetry, thus leading to 
magnetically ordered ground states. Since the 
pairing in ([9]) is a pairing of different bands c CT , 
f a , it docs not have an obvious connection with 
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the spin density operator. However, it can be 
related to local magnetization in a more sophis- 
ticated manner. To illustrate this, we follow 
Ref. [H, and define a local spin operator by 



S(r) 



al(r)cr ss ,a s (r), 



where cr is (a x ,cr y ,a z ), ando^(r) (a ff (r)) repre- 
sents local field creation and annihilation opera- 
tors. An explicit expression for these operators 
is given in Eq. (|A2[) . The local magnetization 
can then be expressed in terms of the spin op- 
erator, 



M(r) = -^£<S(r)) 



(14) 



where the average is taken with respect to the 
dominant ground state obtained from the RG. 
Here, g is the 3-factor, [Ib is the Bohr magne- 
ton, and V is the volume of the system. If we 
expand the local field operator into Bloch waves 
12[) . we obtain 

cj 

E 

ss' A,B K 2 Ki 



xe 



(15) 



where the Bloch wave function is (pa,k(j) = 
e~ lKr \\UA,K (r), r = rj_ + rii, r_L is the out- 
of-plane vector, rn is the in-plane vector, and 
ua,k( t ) is a periodic function with — > r»|| + 
mai+na2, where m, n are integers. 

Let us also introduce a SDW gap function, 



sdw 



3i + Ui 



= E S Z ^ a sA<®° Z ss>f a >sM). (16) 



aa' ;k 



If we confine the system to 2D, setting 
ua,k(t) = S(r±), we can reduce (fT5")) to a sim- 
pler form, 

M(r) = - 2 y^'+^ cexKQ-riiMrx). (17) 

Using this formulation, pairing in the SDW 

channel(cj Cs (fc)(Tj S ,/K's'(fc)) can be easily iden- 
tified by this observable with ordering vector Q. 



Analogously, the FM and FM' pairing chan- 
nels can be identified. Similary, for CDW, we 
introduce the local charge density operator, 

p(r) = ^5>t(r)Mr)>. (18) 

a 

From the mean field Hamiltonian (see Ap- 
pendix [C]), the trial ground state solutions for 
EI and EI' are equivalent to Ref. [29|. In these 
excitonic insulating states, the electrons from 
the conduction band and the holes from the va- 
lence band form bound states. 

Furthermore, FM' and EI' break valley sym- 
metry, i.e. time reversal symmetry, because in 
these ground states, the symmetry exchanging 
K and K' is absent. This can lead to non- 
trivial insulating states^Sr— . Since the mean 
field Hamiltonian in Appendix [Cl does not have 
a clear 'inverted' band gap, we do not conclude 
that these are quantum spin Hall or quantum 
anomalous Hall insulator states. 



B. 



Case II: g , g u g 2 / 



For no n- vanishing values of 50,1,2, much of 
the discussion is similar to the previous section. 
However, since (70.1,2 connect different channels 
in the flow rate equations, they do not give 
simple analytical results that show how the Fj 
evolve. Instead, in this more general case the 
flow rate equations in ([7]) need to be solved nu- 
merically. 

Due to the large parameter space spanned by 
the possible sets of bare couplings, it is impos- 
sible to explore the entire phase diagram. In 
this section, we only select several regions to 
scan, illustrating how finite values of 30,1,2 af- 
fect the results from the previous section. Using 
the bare values from Fig. 21 we scan 30 and g%. 
When 30,1.2 is small, we obtain results very sim- 
ilar to the 30,1.2 = case, with FM occupying 
large regions of the phase diagram. However, 
when 31 becomes large, we instead obtain the 
more complicated phase diagram shown in Fig. 

El 

To form an excitonic insulator, EI and EI' 
order intricately compete with other instabil- 
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FIG. 5. Phase Diagram for a representative choice 
of bare couplings ho = uo = vo = 1, hi = vi = 
Mi = 0.8/io, and /12 = U2 = V2 — gi — O.l/io. The 
phase diagram is determined by monitoring which 
channels divergence first during the RG flow. 



FIG. 6. The flow of go as the bare value of gi varies. 
We set <7o = 32 = O.l/io, and the remaining bare 
couplings are same as Fig. [4] 



itics. Without the scattering processes 50,1,2, 
fine-tuning bare couplings to enhance these in- 
stabilities and suppress the others is inevitable. 
However, introducing finite go.1.2, the flow of 
.90.1,2 significantly affects this result, which can 
automatically enhance or suppress the orders. 
For instance, when g\ is small, the go starts 
flowing towards more positive values (see Fig. 

Consequently, the FM ordering tendency 
is enhanced (which enhances the divergence of 
Tfm)- On the other hand, in Fig. [51 when 
gi becomes large, go flows towards increasingly 
negative values, this suppressing FM order. 
Because of this suppression, EI, EI 1 and SDW 
emerge in the large g\ region as shown in Fig. 

Furthermore, charge density wave order is 
very unlikely to dominate, since Tcdw(0) = 
^SD\v(fi) — 2(v2 +52), and <?2 always grows into 
the positive regime, as long as only repulsive in- 
teractions are considered. We observe that the 
divergence of spin density wave order is always 
stronger than charge density wave order. For 
the same reason, Tfm(0) = TpM'(0) — 2(it 2 + 
g%), and thus FM order is more favorable than 
FM 1 . 

Similarly, superconducting order is not ex- 
pected to dominate for small bare values hi and 
<72- In order to produce dominant BCS insta- 



bilities one needs that h% + 52 > hi + g\ or 
h 2 + .91 > hi + g 2 , such that T SC (0) or T SC '(0) 
is positive. 

To decide which one is the correct ground 
state for BLG is beyond the analysis of this 
paper, because reliable bare coupling constants 
are in general difficult to obtain. The ultimate 
answer will require more experimental input. 



V. SUMMARY 

Summarizing this work, BLG can been mod- 
eled by nearest-neighbor hopping model on a 
bilayer honeycomb structure with 73 ~ 0. A 
general form of interactions between electrons 
can be accounted for without including long- 
range Coulomb interactions. Due to the triv- 
ial topology of the Fermi surface, the RG tree 
level analysis eliminates irrelevant channels and 
greatly simplifies the interacting terms in this 
model. The RG flow rate equation can be calcu- 
lated up to one-loop level in the weak coupling 
limit. 

Instabilities are inevitable if the inter-band 
and inter- valley interactions are nonzero. We 
have investigated each instability, and related 
them to ordered ground states. Specifically, 
we have found competing ferromagnetic (FM, 
FM') , spin density wave (SDW), excitonic in- 
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sulator (EI, EI 1 ), charge density wave (CDW), 
and superconducting (SC, SC) ground states 
in this model. Except SC and SC, all the 
ground states are insulating. Furthermore, val- 
ley symmetry breaking is found in the FM' and 
EI' channels. 

Since the free system with quadratic disper- 
sion is not stable, any small perturbation can 
drive this point toward divergence. Due to 
particle-hole symmetry and perfect nesting of 
the two Fermi points K, K', excitonic insta- 
bilities are expected to arise, which has been 
previously pointed out in Ref. HI and [33l 



VI. DISCUSSION 

Projecting out the orbital field operators 
a\Ka and a-^Ka m © an d taking the spatial 
continuum limit, the non- interacting low energy 
effective model of BLG can be approximated 
by a massive chiral fcrmion modcl£. By sym- 
metry arguments, all possible two-body interac- 
tions can be obtained.— ~— Under this limit, this 
model exhibit a very rich and exotic low energy 
phenomenology because of the newly emerging 
valley and pseudospin degree of freedom. 

In this paper, we utilize the band represen- 
tation point of view which are not necessary 
to impose continuum limit. Projecting out the 
gapped bands, BLG is effectively viewed as a 
conventional two-band model and all its inter- 
actions can be immediately obtained according 
to the band index. In this approach, the in- 
stabilities of BLG is clearly interpreted as the 
peculiar nature of the FS and the low energy 
physics exhibit rich excitonic orders. 

In addition, since we are performing RG 
transformation in band representation, the 
'which-layer — or pseudospin^ symmetry break- 
ing is not obvious, and any 'which-layer' order 
cannot be straight-forwardly observed from the 
order parameters. As discussed in the previ- 
ous section, the order parameters are pairings 
between the c and / bands, thus giving rise 
to a complicated relation for local magnetiza- 
tion and charge density. To extract layer order, 
it would be necessary to know the appropriate 



Bloch wave function or its Wannier represen- 
tation. Comparing to the approach from Ref. 
HcWl4l this is one of the drawbacks of our ap- 
proach. 

Furthermore, the ground states in this pa- 
per have been classified according to their band 
index, and the pseudospin index is implicitly 
contained in ([3]). Therefore, the pseudospin 
symmetry breaking is not made explicit in our 
model. This leads to a different physical inter- 
pretation of the ordered states. Because of this 
reason, not all the of the possible competing 
ground states found in Ref. [lj can be obtained 
by in our study. 

A recent functional renormalization group 
(fRG) study^l has demonstrated the advantage 
of retaining all the lattice structure, and in- 
tegrating out energy modes without ambigui- 
ties. Furthermore, their approach takes into ac- 
count the complication of angular dependence 
in the interactions. Their study has shown an 
interesting "thrce-sublattice CDW instability". 
This instability is quickly disappears as the on- 
site interaction becomes dominant. Since their 
model Hamiltonian is different from the one 
studied in this paper, a direct comparison with 
their results is not straightforward. One of the 
big discrepancies is the predominant FM in- 
stability observed in our approach. This may 
arise because exchange interactions^ were not 
explicitly considered in Ref. [H. 

The results in this paper are valid only of 
the one-loop level. Typically, higher-loop con- 
tributions can be neglected by invoking l/N 
arguments^ However, this type of argument 
is not very strong for this model, because the 
Fermi surface contains just two points, result- 
ing in N — 4 only. Also, the results pre- 
sented here only apply for the weak coupling 
limit. Any strong enough coupling to break 
down the perturbative expansion will invalidate 
the preceding discussion. In the strong coupling 
limit, results from the tree level analysis cannot 
be trusted, and using the same effective action 
as in (|B3[) will not guarantee correct results. 
When the order of the tree and loop diagrams 
is comparable, new effective models and non- 
perturbative approaches may be needed. 
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In parameter space, the point with non-zero 
trigonal warping, i.e. the free part of the action 
is non-analytic. Due to this reason, perturba- 
tivc RG may not work properly. In particular, 
the scaling rule of the fcrmonic field cannot be 
defined, and one does not know whether this 
point is a Gaussian fixed-point or not. There- 
fore, we enforce using 73 = as the fixed-point 
to define the scaling rule, and always treat the 
trigonal warping term as a quadratic perturba- 
tion. 

Furthermore, the instabilities in this paper 
are driven by perfect particle-hole nesting. The 
presence of disorder can destroy this symmetry 
and thus change the phase diagram. Doping 
away from the charge neutrality point, the FS 
becomes a line rather than a few points. In 
this case, the scaling properties of the fcrmionic 
field are different, and the conclusions of our ap- 
proach are no longer valid. As show in a recent 
study of doped monolayer graphene^ func- 
tional RG is a promising alternative method to 
study the BLG doping problem. In the doped 
case, functional RG may superior to our ap- 
proach, since the shape of FS evolves nontriv- 
ially under doping. 

The results of our model are consistent 
with recent current transport spectroscopic 
expcrimentS i 17 i 20 Specifically, a magnetic field 
dependent gap is expected in a ground state 
with excitonic order^i Therefore, magnetically 
order ground states are not a necessary condi- 
tion to exhibit this property. To make connec- 
tions with experiments, the physical properties 
of the ground states discussed in this paper need 
to be analyzed, in particular how these ground 
states respond to external perturbations, espe- 
cially to currents. Also in the experiments, con- 
sidering the effects of disorder and boundaries 
is important. 

We would like to thank Tameem Albash, 
Rahul Nandkishore, Ronny Thomalc, Hubert 
Saleur, Oscar Vafek, and Lorenzo Campos 
Vcnuti for useful discussions, and acknowledge 
financial support by the Department of Energy 
under grant DE-FG02-05ER46240. 



Appendix A: Coupling Constants and the 
Interacting Hamiltonian 

In this appendix, we show how the interacting 
Hamiltonian for the BLG model is constructed. 
To accomplish this, we approximate the Bloch 
wave function by using the 7r z orbital 4>(r) = 
y^J^ze-^ with £ = 1.72oq x ^ 




(Al) 



where N uc is the number of unit cells, eto is the 
Bohr radius, I = c, f bands, R mn = mai + na.2 
denote the center of the unit cell, and t m are 
the basis in the unit cell. K is the in-plane 
crystal momentum of BLG in the first Bril- 
louin zone, i = 1 represents the site 62 with 
ti = (0,0,— c), i = 2 represents the site 02 with 
t2 = (2a, 0, — c), i = 3 represents the site d\ 
with t3 = (2a, 0,0), and i = 4 represents the 
site b\ with t4 = (a, 0,0) (see Fig. [T]). a is the 
lattice spacing 1.57A, and c is the layer separa- 
tion 3.35A. 

Graphenc can be considered as a 2D ma- 
terial, but the electrons still live in 3D real 
space. In order to obtain correct interaction 
terms between electrons, we start from the orig- 
inal Hamiltonian^ (Born-Oppenheimer approx- 
imation is used), which describes the electrons 
with Coulomb interaction in real space repre- 
sentation, 



H fuU = / d 3 rat(r)[— V 2 + V ext (r)]a a (r)+ 
i J d 3 rd 3 r'at(r)a^(r')^ mt (r-r')a ff '(r')a CT (r), 



where aj.(r) (a CT (r)) arc local field operator 
which create (annihilate) an electron at r. V ex t 
is the potential produced by the ions. Vi„t is 
the Coulomb potential between electrons at r 
and r'. 

Now we approximate (the gapped bands are 
neglected) the full operator a a (r) by expanding 
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it into Bloch waves from (|A1[) . Then we have 
at(r)=,^[< K (r)4 x + ^,K(r)4, 



7 if J' 



(A2) 



By using 

r -ri 2 v 2 



2m 



+ V r ext (r)]^ C!/; _ R -(r) = e c j(K)tp c j. K (T) 



and substituting the above equation into Hf u u, 
one can easily obtain i^o in and -ffmt in (0 
and (|A4|) . The coupling constant is determined 

by 



U(K 3 K±K 2 K X 



d 3 rd 3 r'V mt (r - r') 



(A3) 



x <Pi a K a ( r M 4 /f 4 i 2 k 2 {r')tpi lKl (r). 



By substituting (|A1[) into (|A3[) . one can also 
verify that the valley and particle-hole symme- 
tries still hold, yielding only six independent 
coupling constants. Note that projcctinging out 
the gapped bands introduces a hard cutoff, fur- 
ther modifying U, which should behave like the 
original long-range Coulomb interaction^. 

I 



1. Interactions in the BLG Hamiltonian 

From Eqs. ()A3|) and (|5j) we find ten inequiv- 
alent interaction terms, which are not ruled out 
by symmetry and conservation laws,— 

U4(K 3 K 4 K 2 K 1 )fl K3 cl IKi c al K 2 CaK 1 

U^K^Kt)^ ft, K Ja'K 2 fcrKi 
K 3 c a' K t c a'K 2 faK 1 

U 6 (K z K i K 2 Ki)fl K jl, K j a .K*c a K 1 (A4) 

Of these, W 4 and U$ are irrelevant under RG 
tree level, because they vanish at the FS. In the 
following, we will prove W 4 (K, K, K, K) = 0, 
and the proof for different valley combination 
and U$ is similar. 

Using (|A3j) . we have 



Z4(K,K,K,K) = J d 3 xd 3 xV„ K (x)x 
<^* K (x')^ nt (x - x')^ cK (x)< ( 9 c k(x') 



(A5) 



Note that C V K = ^=(1,0,0,1) and C C K = 
^=(1,0,0,-1) at K = K,K'. Now we use Eq. 
(jAll) to write out the Bloch wave function ex- 
plicitly, 



M 4 (K,K,K,K), a ^ E E E |rf 3 xrf 3 x'e !K ' (R "'- 1+R -l»r R — - R ^»^x 

mini m^Bj m 2 n 2 rn' 2 n' 2 

Vi n t(x - x')<L\(x - R mini )$ c (x' - R m ' i „' i )$ c (x' - R TO ^ n ^)$ c (x - R m2 „ 2 ), 



(A6) 



where $„(x) = ^[</>(x - ti) - tp(x - t 4 )], 
$ c (x) = ^=[0(x-ti) + 0(x - t 4 )], and a is 



r 



some constant. Applying changes of variables, 
x —7 x + (ti + t 4 ) and x' — >• x' + (tx + t 4 ), and 
focusing on the <L\(x - R mini )$ c (x - R m2 „ 2 ) 
product term in the integrand. 



$„(x - R mini )$ c (x- R TO2 „ 2 ) 

= j (^( X _ ^-mini - ti) - 0(x - R,„ 2 „ 2 - t 4 ))(0(x - R miTCl - ti) + </>(x - R„ l2 „ 2 - t 4 )) 



-(0(x - R mini + t 4 ) - 0(x - R„ l2 „ 2 + ti))((/)(x - Rmim + 1 4 ) + 0(x - Rm 2 „ 2 + ti)). 
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Next, we perform changes of variables for m and n, m — > M — m and n — > N — n. Therefore, 
R mn — > Rmn ~ Rm.n and this does not affect anything but, 

$„(x — R mi)11 )<I> c (x — R„ l2 „ 2 ) — >-((j>(x — Rmn + Rm ini + ^a) — </>( x — Rmn + Rm 2 n 2 + ti)) 

x((/)(x - Rmn + R mmi + 0(x RmN ~\~ Rm 2 n 2 H" tl)) 



Rmn can be removed by performing changes of 
variables x — > — x+Rmn andx' — > — x'+Rmn- 
Notice that the tt z orbital satisfies 4>(— x) = 
— ^(x). Therefore, we obtain exactly the same 
expression as in (|A6[) . except a minus sign. 
This means Ui must be vanish at FS. Similarly, 
U5 = for the same reason. 

From the result of the RG tree level, all cou- 
plings with small k dependence are irrelevant. 
The first leading non- vanishing term in U^, U§ 
is (D(k). Therefore, the interactions in (|A4[) are 
irrelevant. 



2. Coupling Constant Expansion 

Here we show how to expand the coupling 
constants around K and K'. First we perform a 
unitary transformation by changing from Bloch 
wave to Wannicr representation. The Wannicr 



I 

function is defined as 

B Z - e -iK-R mn 

Wi, mn (r) = — nrr- V>ik{t), (A8) 

K ViVuc 

and 

m,n e iK-n mn 

<Pik(t) = J= Wj,mn(r), (A9) 

where N uc is the total number of unit cells. Also 
from (|A8|) and (|A9|) . we can derive the identity 

J- jK^n = (2k) 2 N uc P (K) (A10) 
111 n 

Note that, e iK - Rm ™ = e i - ff - R """+ G ' Rm " , where 
G is a reciprocal lattice vector. Therefore 8(K) 
in (|A10[) is not exactly the Dirac delta function, 
but equal to a delta function up to a reciprocal 
lattice vector. 

Now, using (|A10[) and (|A3|) . we obtain 



U{K z K i K 2 K 1 ) =^r E E E E e -« 1 .R mWe -ur.-H B1B4e «r r H BJ ^ e «r 1 .R Bini 

uc 7137713 7147714 71 2 7712 Tiimi 

J d 3 xd 3 x'w* l3 , TO3 „ 3 (x)^ 4>m4 „ 4 (x')K ra t( x ~ x')w7 2 , m2 „ 2 (x')w/ limini (x) 



U{K Z K A K 2 K 1 ) =jj2~ E E E Je^W^x 

nC 7237713 7747714 7127712 Til 7711 

e -»3-(B™ 3 „ 3 -R m4 „ 4 ) e iif 2 -(R W2 » a -R m4fl4 ) e tifr(R mi n 1 -R B4 „ 4 ) x 

y d 3 xd 3 x'w* h , TO3n3 (x)u;j 4jro4 „ 4 (x')Vmt(x - x / )u)j a>man2 (x')u;j 1>mini (x) 
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by changes of variables x — > x + R m4fl4 , x' — > x' + R m4 „ 4 and using (|A10[) 
U(K i K 3 K 2 K 1 ) =^-5 2 {K 4 + K 3 - K 2 - Ky) ^ ]T ]T 



^-5 J (K A + K 3 -K 2 -K x ) > > > e 3 m 3- 3e 2 m ^ e 

n'^ra'^ n' 2 m' 2 n^m^ 

d 3 xd 3 x'w* l3fi0 (x)w* Ii m , n , (x')^„t(x-x')u; /2im ^ i (x')w /ljmi „ ; (x) 

(All) 



The coupling constant expansion can 
be achieved by expanding e i 3 Rm 3"3 , 

iK 2 -R m , i iKi-~R.ii . „ /TTTTT 

e m 2"2 ; e ii in -hq. (jAllj). 

S(K 4 + K 3 — K 2 — Ki) means that mo- 
mentum is conserved up to a reciprocal lattice 
vector. This allows Umklapp processes. Mo- 
mentum conservation emerges because of the 
in-plane lattice translational invariance in the 
system. 

Appendix B: Renormalization Group and 
Tree Level Analysis 

This section summarizes the details needed 
for the pcrturbative renormalization group 
analysis. 



1. Action of the Model Hamiltonian 



The RG transformation is performed using 
the path integral formalism. Therefore, the 
model Hamiltonian from Section 2 should be 
rewritten into action form. S = J drC = 
So + Si n t, where So is the free action and Smt 
contains the interaction terms. The deriva- 
tion is tedious, and one needs to introduce co- 
herent states of the creation and annihilation 
field operators^. However, the result is sim- 
ple, which can be achieved by replacing the field 
creation and annihilation operator by Grass- 
mann fields. Namely, c) a — > c a — > ip a and 
ft -* Xa, f -> X<r- Therefore, the S and S lnt 
can be written as 



duj f d 2 K 



S ° = 2~ \k-k\<a 7^MK,u)(iu ^ e c {K))i' a {K,u) + XaiK.uj)^ - ^{K^XaiK.uj) 

J-oo * • / |^ K <|< A ^ n ) 



Si: 



l r A [°° dui f d 2 K 

i=l J -°° J \K z -K'\<A y ' 



+ U 1 (K 3 K 4 K 2 K 1 )^(K 3 ,U3)^AKi^4)xAKa,(^)Xa(K 1 ,u 1 ) ( B1 ) 
+ U 2 {K 3 K4,K 2 K-{)-4) a {K 3 , uj 3 ) X *> {Ki,ui)x*> i^ 2 , u 2 )M K i, wi) 
+ U i (K 3 K 4 K 2 Ki)xa(K 3 ,uj 3 )^ crl (K 4 ,u} 4 )xa'(K 2 ,uj 2 )^ cr (Ki,u}i)} 
+ [exchange (tp -H> x)] 



Now we can perform the RG analysis. First, This choice will determine the scaling proper- 
Sq is chosen to be the fixed point in the theory. 
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ties of u! and if), which will be discussed in the 
following section. 



Scaling Properties and Effective Action 
at the Tree Level 



Since we interested in low energy limit, only 
the energy modes in the vicinity of Fermi points 
are considered (see Fig. 



2(a) |. Expanding 



e c j(K) around K and K', and combining with 
the results from @, 



E 

q=K,K' 



00 du 

-or, 27T 



fc|<A 



d 2 k 



i> aa (k, - ^-k 2 )ip ar7 (k,uj) 

2 



(B2) 



We introduce a short hand notation 
ip aa (k,u)) = i; a (a + fc,w), tp a<J (k,uj) = 
4><j{a + k,u>), Xaa(k,uj) = Xa(a + k,oj), 
Xaa(k,uj) = x<r( a + k,ui), and a = K,K' is 
known as the 'valley' degree of freedom. Valley 
index is similar to L (left) and R (right) index 
in the one dimensional case^ 

The RG transformation is simply integrat- 
ing out the high energy modes of ip aa {k,ui), 
ip aa (k,u}) and x Q(T (fc,w), Xaa(k,uj) which lie 



within the thin shell, A, in Figure 2(a)| and 
considering how these modes affect the low en- 
ergy theory. After integrating out, only \k'\ < 
A— dA = A/s modes remain in the theory. In or- 
der to evaluate what has changed from the orig- 
inal theory, k' must be rescaled (k' = sk) back 
to the original phase space such that |fc| < A. 
Since So is the fixed point, this requires that uj, 
ip a:Cr (k,Lu) and Xa,a{k,oS) must be rescaled, 

ui' = s 2 ui, 
1 p' aa (k',u; , ) = s- 3 ^ cta (k'/s,uj), 
X' aa (k',u') = s-^ Xa a{k' /s,uj). 

With this scaling relation, one can now ask 
how the coupling constants Uq, Hi, U2, U3, 
U/± and U$ scale under the RG transformation. 
Again, we use Eq. (|A11[) to expand couplings 
around K and K'. By enforcing momentum 
conservation, only the constant term in the ex- 
pansion do not renormalize to zero (marginal 
under tree level). 

Note that Si„ t remains unchanged when 
K4 -H* A3 and K2 <-> A'i simultaneously. In 
addition, using time reversal symmetry (val- 
ley symmetry) in the model, hence exchanging 
K <-> K' in (Table [IJ will not produce another 
set of independent coupling constants. Thus we 
obtain Si n t at the tree level, 
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■5 [ft 

i=l 



|fc»|<A 



(2tt) 2 - 



(27r) 2 (5 2 (fci + fa — fa — /c 4 )27T(5(a;i + cj 2 — w 3 — cj 4 ) x { 



Wk<t (&3 j W 3 ) ^Kcr' (&4 , W 4 ) 1pKa> (fa , W 2 ) V'Kcr (fcl , U>i) 

+/ii^Kct(*3, w 3 )i/)kv'(^4, w 4 )V'K'o-'( fc 2, w 2 )^Kff(fci, ^i) +exchange (K -h- K') 

+h2$K>a (fa , W 3 )-0Ka' (&4 1 ^KV (fe i ^2 )fc (fcl , 
SoV'Ko- (fe , W 3 )V>K<7' (&4i W 4 )XK(T' (&2, W2)XKa (&1 , W 4 ) 

+ffi?/'K CT (A: 3 ,w 3 )i/5K'(T'(fc4,w 4 )xK' C r'(fc2, oj 2 )XK<r(fa, u>i) +exchange (K o K') 

+92lpK.'a (fa , W 3 )^ K <t' (fa , W 4 )XK'<T' (&2 , ^2 )XK<x (fa , Wl ) 
MO V'Kcr (fc 3 , W 3 )xk<t' (&4 , W 4 )xk<t' (fa , W 2 )l/>Ka (fa , Wl) 

+«i^Kafe,w 3 )xKv(fc4,W4)XKv(fc2, ^^KaOi, ^i) +exchange (K «-» K') 

_+W2^KV (fa , W 3 )XK<T' (fa , W 4 )XKV (fa , ^2 ) V"Ka (fa , Wl) 
W0XK<t (fa 1 W 3 )^K(r' (fa , W 4 )XK CT ' (fa , U 2 )lpKa (fa , W X ) 

+«iXK C r(fa,w 3 )V> K ' C r'(fa, w 4 )xK'<r'(fa,w 2 )V'K<7(fa,wi) +exchange (K o- K') 

_+W2XK'<r(fa, W 3 )^K<t' (fa, ^4)XK' CT ' (fa, U>2)lpK<7 (fa > ^l) 

+ [exchange O x)] }■ 



(B3) 



Appendix C: Mean Field Analysis of the 
Ground States 

In this section, we summarize the mean 
field analysis for the ground states FM, FM', 
SDW, EI, EI 1 , and CDW. The idea of the 
mean field approximation^ is to guess a trial 
ground state which can minimize the total en- 
ergy of the many-body system. With a given 
trial ground state, the original Hamiltonian can 
be approximated by an effective quadratic mean 
field Hamiltonian which can be solved by self- 
consistent diagonalizing. 

The procedure will be briefly shown in the 
following. First, let the mean field Hamiltonian 
be 

HMF = H o + Hp a i r 

Ho is the free Hamiltonian given in Eq. ((4]), and 

Hpair ^ ^■a<x';a(7 t (fa k )c a(T j i f a , l j / ^'-\-h.C. 

aa' \kk' 

For simplicity, since the coupling constants 
are independent of k and k', we assume 



A^(fafa) ~ A^, iW 5 2 (fc - fa) . Therefore, 

the pairing gap function A^, , (order param- 
eter) is given by 



A 



fin 



Gfr, 



A /m ,' , 



Gfr 



A 



sdw 



G s 



a: 



G e i 

A 68 ' , , 

aa' ;ss' 
G e i' 
^cdw 



G, 



cdw 



k 
k 
k 

^^2( C as,k^ aa ' ® <*s«'/a'«',k) 
k 

J2( C L,k T * a > ® S ss'fa>s',k) 
k 



where Gj = Fj(0) is a linear combination of the 
bare coupling constants listed in Table [TTJ To 
find the trial ground state, a new quasi-particle 
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(mixture of c, v bands) is introduced which di- 
agonalizes the mean field Hamiltonian, 

aa \<7<7 a a \oo p 

Voter, k — ^fc C-occr,k Ja'cr,k 

\ ota. \aa' . aa ;crcr r 

^a<r,k — ^fe C-acr,k T" ^ Ja'cr,k 

Where |u(fc)| 2 + |u(fc)| 2 = 1. 

{A^.A^fc')} = S aotl 5^5 2 {k-k') 

The other commutation relation are zero. Then, 
the diagonalizcd Hamiltonian can be written as 

H MF =J2 E ( k ){vL(k)VaAk)-^L(k)X a Ak) } 



tr.k 



where, 



and 



E(k) 



6(fc)| = 



|A(°)| 2 , 



A 



(o) 

aa' \crcr' 



V2 L 

\ I aa :crcr 



\ E(k)J ' 



e (k) = |e±(fc)| is given by Eq. © 

Therefore, using a Hartree-Fock state as the 
trial ground state for the quasi-particles, we 
have U a ,k ALWI°) = l*AO)>, where |0) is the 
state with no particles (vacuum). To calculate 
, one can apply l^^toj) to calculate the or- 
der parameter and obtain the gap equations, 



A 



(0) 



E 

k 



GA 



(o) 



^)| 2 + |A^ 



(0) 12 
aa 1 I 



OLOL :<7<7 



V2 



<k) 
E(k) 



which are then solved self-consistcntly. 



A. H. Castro Neto, F. Guinea, 
Peres, K. S. Novoselov, and A. 



|Rev. Mod. Phys. 81, 109 (2009)| 
D. E. Sheehy and J. 



N. M. R. 
K. Geim, 

Schmalian, 



Phys. Rev. Lett. 99, 226803 (2007) | 



M. A. H. Vozmediano, Royal Society of London 
Philosophical Transactions Series A 369, 2625 
(2011), |arXiv:1010.5057 [cond-mat.str-el]| 
D. C. Elias, R. V. Gorbachev, A. S. May- 
orov, S. V. Morozov, A. A. Zhukov, P. Blake, 
L. A. Ponomarenko, I. V. Grigorieva, 



Novoselov, F. Guinea , 
Nat Phys 7, 701 (2011)| 



and A. K. 



K. S. 

Geim, 



A^ Cortijo, F. Guinea, and M. A. H. 

Vozmediano, ArXiv e-prints (2011), 
arXiv: 1 1 12 . 2054 [cond-mat .mes-hall] | 
V. N. Kotov, B. Uchoa, V. M. Pereira, F. Guinea, 

and A. H. Castro Neto, ArXiv e- prints (2010), 
arXiv:1012.3484 [cond-mat.str-el]) 

E. McCann and V. I. Fal'ko, 



Phys. Rev. Lett. 96, 086805 ( 2006) 



H. Min G. Borghi, M. Polini, and A. H . Mac- 
Donald, |Phys. Rev. B 77, 041407 (2008) 
R. Nandkishore and L 



Nandkishore 
Phys. Rev. Lett. 104 



Phy; 



156803 (2010) 



Vafek 



and 



K. 



s. Rev. 

F. Zhan g, H. Min, M. Polini, and A. H 
Donald, |Phys. Rev. B 81, 041402 (2010) | 
Y. Lemonik, L L. Aleiner, C. Toke, and V. I. 
Fal'ko, j|Phys. Rev. B 82, 201408 (2010)| 
R. E. Throckmorton and 

O. Vafek, ArXiv e-prints (2011), 



B 81, 041401 (2010) 



Levitov, 
Yang, 
Mac- 



arXiv : 1 1 1 1 . 20 76 [cond-mat . str-el] 



E. 



V. 
and 

|arXiv:1206.0288 [cond-mat.str-el] 



Cvetkovic, 
O. Vafek, 



R. 

ArXiv 



Throckmorton , 



e-prints 



R 



Nandkishore and 
Phys. Rev. B 82, 115124 (2010)| 



(2012), 



Levitov, 



A. S. Mayorov, 
Kruczynski, R. V. 
A. Zhukov, S. V. 



D. C. Elias, 
Gorbachev, T. 
Morozov, M. I. 



M. Mucha- 
Tudorovskiy, 
Katsnelson, 



17 



V. I. Fal'ko, A. K. Geim 
|Science 333, 860 (2011J 



and K. S. Novoselov, 



J. Velasco, L. Jing, W. Bao. Y. Lee, P. Kratz, 
V. Aji, M. Bockrath, C. N. Lau, C. Varma, 
R. Stillwell, D. Smirnov, F. Zhang, J. Jung, and 
A. H. MacDonald, |Nat Nano 7, 156 (2012) 
J. Martin, B. E. Feldman, R. T. 
Weitz, M. T. Allen, and A. Yacoby, 



Phys. Rev. Lett. 105, 256806 (2010) | 



R. T. Weitz, M. T. Allen, B. E. Feldman, J. Mar- 



tin, and A. Yacoby, Scie nce 330, 812 (2010) 
F. Freitag, J. Trbovic, 

M. Weiss, and C. S chonenberger, 

|Phys. Rev. Lett. 108, 076602 (2012")] 
Including trigonal warping term at the beginning 
will make the kinetic term non-analytic. Due to 
this reason. 

J. Nilsson, A. H. Castro Neto, N. M. R. Peres, 



and F. Guinea, [Phys Rev. B 73, 214418 (2006) 
T. Stroucken, J. H. Gronqvist, and S. W. Koch, 



Phys. Rev. B 84, 205445 (2011)| 

R. Shankar, |Rev. Mod. Phys. 66, 129 (1994)| 
A. V. Chubukov, D. V. Efremov, and I. Eremin, 
|Phys. Rev. B 78, 134512 (2008) 



R. Nandkishore^ II S. Levitov, and A. V. 
Chubukov, |Nat Phys 8, 158 (2012)] 
A. Chubukov, Annual Review of Condensed 
Matter Physics 3, 57 (2012). 



M. 



R. Brydon and C. Timm, 



|Phys. Rev. B8 0, 17440 1 (2009) 

D. Jerome, T. M. Rice, and W. Kohn, 



Phys. Rev. 158, 462 (1967) 



Kane and 



E. 



C. 

|Phys. Rev. Lett. 95, 226801 (2005) 



Mele, 



S. Rag hu, X.-L. Qi, C. Honerkamp, and S.-C. 
Zhang, |Phys. Rev. Lett. 100, 156401 (200"8)] 
X.-L. 



Qi 



and 



S.-C. 



Rev. Mod. Phys. 83, 1057 (2011) 



R. 



Nandkishorc 



Rev . B 82, 1 15431 
Scherer 



a nd 
(2010)| 



Phys. 

W. M. Scherer, ST Uebelacker, 

C. H onerkamp, ArXiv e-prints 
|arXiv:1112.5038 [cond-mat-str-el]] 
H. Wang and V. W. 
|Phys. Rev. B 85, 075438 (2012)| 

M. L. Kiesel, C. Piatt, W. 

D. A. Abanin, and R. 



Zhang, 

Levitov, 

and 
(2011), 



Scarola, 

Hanke, 
Thomale, 



Phys. Rev. B 86, 020507 (2012; 



E. W. Fenton, |Phys. Rev. 170, 816 (1968)| 
A. Altalnd and B. D. Simons, Condensed Mat- 
ter Field Theory, 2nd ed. (Cambrige University 
Press, 2010). 

X.-G. Wen, Quantum Field Theory of Many- 
Body System (Oxford University Press, 2007). 



